A Lattice Boltzmann Model for Wave and Fracture phenomena 
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We show that the lattice Boltzmann formalism can be 
used to describe wave propagation in a heterogeneous media, 
as well as solid-body-like systems and fracture propagation. 
Several fundamental properties of real fractures (such as prop- 
agation speed and transition between rough and smooth crack 
surfaces) are well captured by our approach. 
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Lattice Boltzmann (LB) models jj],D are dynamical 
systems, discrete in time and space, aimed at simulating 
the behavior of a real physical system in terms of local 
density of fictitious particles moving and interacting on 
a regular lattice. The density distribution functions are 
denoted fi(f,t) where r refers to the lattice site, t the 
iteration time while the subscript i labels the admissible 
speed of motion Vj (e.g. along the main lattice direc- 
tions). The value i = corresponds to a population of 
rest particles with vq = 0. 

Lattice BGK models || have been used successfully 
to simulate fluid dynamics and complex flows . The 
same approach can be adapted to model wave propaga- 
tion in a heterogeneous media, where propagation speed, 
absorption and reflection can be adjusted locally for each 
lattice sites. We show how such a model can be derived 
and applied to the study of fracture propagation 

Lattice BGK models are characterized by the following 
dynamics 

Mr + riU + r) - Mr, t) = j (jf ) (r, t) - fr(r, <)) 



(1) 



where r is the time step, (r, t) the so-called local equi- 
librium distribution and £ a relaxation time. The func- 
tion f- ' is the key ingredient for it actually contains the 
properties of the physical process under study: this is the 
distribution to which the dynamics spontaneously relaxes 
and which is, therefore, intimately related to the nature 
of the system. 

Wave phenomena, whether mechanical or electromag- 
netic derives from two conserved quantities \& and J, 
together with time reversal invariance and a linear re- 
sponse of the media. The quantity ^ is a scalar field and 
J its associated current. For sound waves, *F and J are 
respectively the density and the momentum variations. 
In electrodynamics, >F is the energy density and J the 
Poynting vector B. 



The idea behind the LB approach is to "generalize" a 
physical process to a discrete space and time universe, 
so that it can be efficiently simulated on a (parallel) 
computer. For waves, this generalization is obtained 
by keeping the essential ingredients of real phenomena, 
namely conservation of \P and J, linearity and time rever- 
sal invariance. Thus, in a discrete space-time universe, a 
generic system leading to wave propagation is obtained 
from cq. (|I|) by an appropriate choice of the local equi- 
librium distribution 



J 
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if i ^ 0, and f£ 0) = a * (2) 



where v is the ratio of the lattice spacing to the time 
step, and \& and J are related to the fiS in the standard 
way: 'F = Xii an d J = Ylifi^i- Note that, here, we 
make no restriction on the sign of the /jS which may well 
be negative in order to represent a wave. 

As opposed to hydrodynamics [EJ , is a linear func- 
tion of the conserved quantities, which ensures the su- 
perposition principle. The parameters a, b and do are 
chosen so that >F = Y^ifi anc ^ J = Y^i^ifi ; > wn ich 
ensures conservation of *F and J. For a two-dimensional 
square lattice (D2Q5 to use the terminology of 0) we 
find cio + 4a — I and 6=1. The freedom on the value 
of ao can be used to adjust locally the wave propagation 
speed. Time reversal invariance is enforced by choos- 
ing £ = 1/2 as can be easily checked from eq. (|I|) with 
J — > —J and * — > * in relation (§). Note that the 
D2Q5 lattice is known for giving anisotropic contribu- 
tions to the hydrodynamic equations. These terms are 
not present in our wave model because they appear with 
a vanishing coefficient when £ = 1/2 (see eq. (f|). 

In hydrodynamic models, £ = 1/2 corresponds to the 
limit of zero viscosity, which is numerically unstable. In 
our case, this instability does not show up provided we 
use an appropriate lattice. Indeed, in the D2Q5 lattice, 
our dynamics is also unitary H which ensures that ff 
is conserved. This extra condition prevents the fiS from 
becoming arbitrarily large (with positive and negative 
signs, since ^ is conserved). This is no longer the case 
with the D2Q9 lattice, where numerical instabilities de- 
velop for our wave dynamics. This observation may shed 
some light on the origin of the numerical instabilities ob- 
served in hydrodynamic models. 

Note that dissipation can be included in our micro- 
dynamics. Using £ > 1/2 allows us to describe waves 
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with viscous-like dissipation. This makes sense with the 
hexagonal lattice D2Q7, where no stability problem oc- 
curs when £ = 1/2 and no anisotropy problem appears 
when the viscosity is non-zero (£ > 1/2). Below we shall 
propose another way to include dissipation in the square 
lattice model, which will be appropriate to our purpose 
of modeling fracture propagation. 

The multiscale Chapman-Enskog expansion || can be 
used to derive the macroscopic behavior of \& and J when 
the lattice spacing and time step go to zero. We obtain 



ft* + dpjp = 



d t J a + 2av 2 d a y 



(2£ - 1) arv 2 d a div J - -^T a p lS dpd 1 J s 



(3) 

= 
(4) 



where T a p 1 s = ^ Vi a VipVi 7 Vis and summation over re- 
peated greek indices (which label the spatial coordi- 
nates) is assumed. With £ = 1/2 equation (||) becomes 
dtJ a + 2av 2 d a 1 $ l = 0. When combined with equ. (||), we 
obtain 

ft 2 * - 2au 2 V 2 * = 

which is a wave equation with propagation speed c = 
v\/2a (note that v is the speed at which information trav- 
els). As mentioned previously, the propagation speed c 
can be adjusted from place to place by choosing the spa- 
tial dependency of a. Provided that ao + 4a = 1 and 
a o > (required for stability reasons), the largest pos- 
sible value is a = 1/4 and corresponds to a maximum 
velocity cq = v/y/2. Therefore media with different re- 
fraction indices n — c /c~ l/(2y/a) can be modeled. 

Perfect reflection on obstacles can be included by mod- 
ifying the microdynamics to be fi (r + TVi ,t + r) = 
—fi>(r,t) on mirror sites, where i' is defined so that 
Vi> = —Vi, i.e. the flux bounces back to where they came 
from with a change of sign. Absorption on non-perfect 
transmitter sites can be obtained by modifying the con- 
servation of <3> to J2i fi — where < fi < 1 is an 
attenuation factor. This modifies a — > fia and ao — » fiaQ. 
Finally, by substituting (|2|) into (|l|) and using the ex- 
pression of a and ao in terms of c, free propagation with 
refraction index n(r), and partial transmission and re- 
flection can be expressed as 



fi(r + TVi, t + r) 



f (?,t + T) = 2/i 



n 2 - 1 



*-/o(r,i) 



(5) 



In this equation, /i = corresponds to perfect reflection, 
fx = 1 to perfect transmission and < /i < 1 describes a 
situation where the wave is partially absorbed. A partic- 
ular version of our LB wave model has been successfully 



validated by the problem of radio wave propagation in a 
city§. 

The idea of expressing wave propagation as a discrete 
formulation of the Huygens principle has been considered 
by several authors ■ Not surprisingly, the resulting 

numerical schemes bear a strong similarity to ours. Nev- 
ertheless the context of these studies is different from ours 
and none have noticed the existing link with the lattice 
BGK approach. Models of refs. [p|,|l0t use a reduced set 
of conserved quantities, which may not be appropriate in 
our case. Other models Jl~T[ ] consider wave propagation in 
a LB approach, but with a significantly more complicated 
microdynamics and a different purpose. 

In what follows, we show how our LB dynamics can 
model a solid body and capture the generic feature of 
crack propagation. Whereas LB methods have been 
largely used to simulate systems of point particles which 
interact locally, modeling a solid body with this ap- 
proach (i.e modeling an object made of many particles 
that maintains its shape and coherence over distances 
much larger than the interparticle spacing) has remained 
mostly unexplored. A successful attempt to model a 
one- dimensional solid as a cellular automata is described 
in |l2j . The crucial ingredient of this model is the fact 
that collective motion is achieved because the "atoms" 
making up the solid vibrate in a coherent way and pro- 
duce an overall displacement. This vibration propagates 
as a wave throughout the solid and reflects at the bound- 
ary. 

A 2D solid-body can be thought of as a square lat- 
tice of particles linked to their nearest neighbors with a 
spring-like interaction. Generalizing the model given in 
p2| requires us to consider this solid as made up of two 
sublattices. We term them black and white, by analogy 
to the checkerboard decomposition. The dynamics con- 
sists in moving the black particles as a function of the 
positions of their white, motionless neighbors, and vice- 
versa, at every other steps. 



Let us denote the location of a black particle by ■ 



(xij,yij). The surrounding white particles will be at 
positions ri-ij, ^i+ij, an d We define the 

separation to the central black particle as (see figure [l]) 



fx(i,j) 


~ r ij 




+ h) 




= r i,j 


- (n,j-i 


+ u) 


/3(£,i) 


= r i-j 


- (n+i,j 


-h) 


U(i,j) 


— r i,j 


- {n,j+i 


-u) 



(0) 



where the fi are now vector quantities, and h = (ro,0) 
and u = (0, ro) can be thought of as representing some 
equilibrium length tq of the horizontal and vertical spring 
connecting adjacent particles. With this formulation, the 
coupling between adjacent particles is not given by the 
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Euclidean distance but is decoupled along each coordi- 
nate axis (however, a deformation along the x-direction 
will propagate along the y-direction and conversely). 
This method makes it possible to work with a square 
lattice, which is usually not taken into account when 
describing deformation in a solid because, with the Eu- 
clidean distance, the y-axis can be tilted by an angle a 
without applying any force. The breaking of the rota- 
tional invariance is expected not to play a role in the 
fracture process we shall consider below. 




FIG. 1. Illustration of the way the fiS are denned. The 
cross indicates the location of the geometrical center of mass 
of the four white particles. At the next iteration, the black 
particle jumps to a symmetrical position with respect to this 
point. 

The motion of all black particles is obtained by up- 
dating the above fas by equ. (||), with n = 1 and for 
i > 0. The local value of $ (which is conserved by the 
dynamics) is then interpreted as the momentum p of the 
corresponding particle. The new location (i+1) of par- 
ticle ij is thus obtained as fij(t+ 1) = ry(t) — (l/my)$, 
where m^- is a mass associated with the particle. Next, 
the quantities / are propagated to the neighbors and in- 
terpreted as the deformations seen by the white particles, 
which then move according to the same procedure. 

A geometrical interpretation of this dynamics is given 
in fig. [l| for n = 1, a bulk particle moves to a symmetric 
location with respect to the center of mass of its neigh- 
bors, (1/4) [Fi-ij + h + f i+ i t j — h + fij-i + u + Fij+i - u] . 
Thus, our model corresponds to looking at particles in 
a harmonic potential, at the particular time where they 
have maximal displacement and no velocity Note also 
that, in this interpretation J x and J y give, respectively, 
the separation between the left-right and up-down neigh- 
bors. 

At the boundary of the domain, or for particles with 
broken bonds, a different rule of motion has to be con- 
sidered. The above interpretation, in terms of a jump 
with respect to the center of mass of the neighbors, gives 
a natural way to formulate the dynamics when less than 
four links are present. 

Our original wave paradigm includes five fields. In the 
case of a solid body, the fifth quantity /o can be added in 
the model as an internal degree of freedom. This is useful 
in order to describe solids with different sound speeds 
c and to then see how the fracture propagation speed 



relates to c. The particle motion is still determined by the 
local momentum \& but, now, there is no more distinction 
between white and black particles. 

Our goal is now to show that our LB wave model can 
be used to describe a fracture process. Fracture is a phe- 
nomena for which no definite theory is available |l3| ] and 
a simple model is certainly useful to help understanding 
generic properties. 

At the level of our description, a fracture is easily intro- 
duced. A link may break locally if its deformation is too 
large. Here we consider the energy stored in a link as the 
quantity determining the breaking. The energy Ej. (r , t) 
of each link k, (k = 1 to 4) at site r and time t is defined 
as Ek = (l/4)/o + fk- Since, as mentioned above, our 
dynamics is unitary, the total energy E to t — J2k fEk(r) 
is conserved until a link breaks. Note that, according 
to our interpretation, the microscopic energy is only of 
potential type at integer time steps. 

The breaking rule we impose is as follows: a link k 
breaks if the corresponding Ek is larger than a given 
threshold e(r) which may, in principle, depend on the 
position (local defects). Particles with one or more bro- 
ken links then behave like particles at the boundary. 

The total energy E tot can be written as the kinetic 
energy of the center of mass Ekin, phis an internal energy 
E- lnt = E to t — -Ekin- The kinetic part is computed as 
Ekin = (l/2M)[^pr/] 2 , where M is the total mass. 

The second contribution, E m t, can be set proportional 
to a temperature T, using the equipartition theorem. 
In the initial configuration, T is typically introduced by 
adding a noise of standard deviation y/T to the rest po- 
sition of each atom. 











v 



FIG. 2. Contour plot of the energy E in a 100 x 100 square 
solid under a preloaded (mode I) stress (S = 2.2 x 10 -2 ). The 
contour line mark the levels e/E — 10 x £, with I = 1, 10. 

A typical experiment which is performed when study- 
ing fracture formation is to apply a stress by pulling in 
opposite directions the left and right extremities of a solid 
sample. To achieve a static stress, the solid is prepared in 
a configuration where the re-spacing between the atoms is 
increased to the value (l+S)ro where S is called the stress 
factor. Both left and right extremities are not allowed to 
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move. The initial temperature may be different from 
and a small notch (artificially broken links) is made in 
the middle of the sample to favor the apparition of the 
fracture at this position. In the fracture community, this 
is known as mode I loading JT3] . 

Figure [2] shows the stationary spatial distribution of 
stress around the notch given by the energy contour lines 
and obtained from our simulation in a case where the 
fracture does not propagate across the sample but stops 
after a few steps. In figure || we can see the result of 
two fracture experiments where the crack spontaneously 
propagates through the sample after being initiated ar- 
tificially. The dissipation coefficient /1 in (^|) turns out 
to be essential and eventually distinguishes the two cases 
presented here. Attenuation prevents the reflection of 
too much energy from the boundaries toward the crack. 
With /j, — 0.91 the dissipation of energy is high enough 
to limit the acceleration of the crack below some criti- 
cal speed and the crack remains smooth. On the other 
hand, with /i = 0.96 the crack accelerates above this crit- 
ical speed and instabilities appear: the crack progresses 
while making micro-branching. 



were performed on samples with c = l/v2, but we ob- 
serve a similar behavior using another values of c. Thus, 
a critical fracture speed (l3) of about c/2 for the smooth- 
branching transition is well captured in our model. This 
non-trivial result is promising since no simple statistical 
models are yet available to describe a fracture process. 
The fact that our dynamics is based on a description at 
the particle level makes the connection with real experi- 
ment possible and leaves much flexibility to adjust locally 
some parameters. Yet, our model is significantly simpler 
than a molecular dynamics approach. As quantitative 
experimental results are now possible at a microscopic 
scale (e.g. scaling in the apparition of micro-cracks be- 
fore the main crack |p"5|| ), our model could contribute to 
highlight the essential mechanisms of fracture. Another 
interesting check would be to measure, both numerically 
and experimentaly the relation between roughness and 
fracture speed. 

We acknowledge support from the Swiss National Sci- 
ence Foundation. 
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FIG. 3. Fracture tip location as a function of time (itera- 
tions) for two types of cracks; for comparison, straight lines 
corresponding to half the speed of sound are drawn. The 
fastest crack corresponds to a branching situation (see fig- 
ure §), obtained for S = 0.03, y, = 0.96, e = 0.0058. The 
slowest crack is smooth (see figure ^) and is produced with 
S = 0.03, fi = 0.91, e = 0.0058. 




FIG. 4. Map of the broken links for the smooth (left) and 
branching cracks (right) described in figure |j| 

As can be seen in figure |^, the limiting speed is around 
half of the speed of sound in the sample. The simulations 
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